%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cs\_turbulence\_ke}\label{ap:turbke}

\hypertarget{cs\_turbulence\_ke}{}

\vspace{1cm}
%-------------------------------------------------------------------------------
\section*{Function}
%-------------------------------------------------------------------------------
The purpose of this function is to solve the system of equations of
$k$ and $\varepsilon$ in a semi-coupled manner.\\
The system of equations solved is the following:

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho\frac{\partial k}{\partial t} +
\dive\left[\rho \vect{u}\,k-(\mu+\frac{\mu_t}{\sigma_k})\grad{k}\right] =
\mathcal{P}+\mathcal{G}-\rho\varepsilon+k\dive(\rho\vect{u})
+\Gamma(k_i-k)\\
\multicolumn{1}{c}{+\alpha_k k +\beta_k}\\
\displaystyle
\rho\frac{\partial \varepsilon}{\partial t} +
\dive\left[\rho \vect{u}\,\varepsilon-
(\mu+\frac{\mu_t}{\sigma_\varepsilon})\grad{\varepsilon}\right] =
C_{\varepsilon_1}\frac{\varepsilon}{k}\left[\mathcal{P}
+(1-C_{\varepsilon_3})\mathcal{G}\right]
-\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k}
+\varepsilon\dive(\rho\vect{u})\\
\multicolumn{1}{c}{+\Gamma(\varepsilon_i-\varepsilon)
+\alpha_\varepsilon \varepsilon +\beta_\varepsilon}
\end{array}\right.
\end{equation}

$\mathcal{P}$ is the term of production by mean shear stress:
\begin{displaymath}
\begin{array}{rcl}
\mathcal{P} & = & \displaystyle -\rho R_{ij}\frac{\partial u_i}{\partial x_j}
= -\left[-\mu_t \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)
+\frac{2}{3}\mu_t\frac{\partial u_k}{\partial x_k}\delta_{ij}
+\frac{2}{3}\rho k\delta_{ij}\right]
\frac{\partial u_i}{\partial x_j}\\
& = & \displaystyle \mu_t \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}
-\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})\\
& = & \displaystyle \mu_t \left[
2\left(\frac{\partial u}{\partial x}\right)^2+
2\left(\frac{\partial v}{\partial y}\right)^2+
2\left(\frac{\partial w}{\partial z}\right)^2+
\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2+
\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)^2+
\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2
\right]\\
\multicolumn{3}{r}%
{\displaystyle -\frac{2}{3}\mu_t(\dive\vect{u})^2-\frac{2}{3}\rho k \dive(\vect{u})}
\end{array}
\end{displaymath}

$\mathcal{G}$ is the gravity production term:
$\displaystyle
\mathcal{G}=-\frac{1}{\rho}\frac{\mu_t}{\sigma_t}
\frac{\partial\rho}{\partial x_i}g_i$

The turbulent viscosity is
$\displaystyle \mu_t=\rho C_\mu\frac{k^2}{\varepsilon}$.

The constants are:\\
$C_\mu=0,09$ ;
$C_{\varepsilon_2}=1,92$ ; $\sigma_k=1$ ; $\sigma_\varepsilon=1,3$\\
$C_{\varepsilon_3}=0$ si $\mathcal{G}\geqslant0$ (stratification unstable) et
$C_{\varepsilon_3}=1$ si $\mathcal{G}\leqslant0$ (stratification stable).

$\Gamma$ is a possible mass source term (such that the mass conservation equation becomes
$\displaystyle \frac{\partial \rho}{\partial t}+\dive(\rho\vect{u})=\Gamma$).
$\varphi_i$ ($\varphi=k$ or $\varepsilon$) is the value of $\varphi$
associated with the injected or extracted mass. In the case where we remove
mass ($\Gamma<0$), we necessarily have $\varphi_i=\varphi$.
Similarly, when we inject mass, we also often specifies $\varphi_i=\varphi$.y
In these two cases, the term disappears from the equation. In the following sections,
we will qualify as {\em forced injection} the cases where we have $\Gamma>0$ and
$\varphi_i\ne\varphi$.

$\alpha_k$, $\beta_k$, $\alpha_\varepsilon$, $\beta_\varepsilon$ are possible
user source terms, leading to partial implication,
imposed if necessary by the function \fort{cs\_user\_source\_terms}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discretization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The resolution is done in three steps, in order to partially couple the two variables
$k$ and $\varepsilon$. For simplicty, let us rewrite the system as follows:

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho\frac{\partial k}{\partial t} =
D(k) + S_k(k,\varepsilon)+k\dive(\rho\vect{u})+\Gamma(k_i-k)+\alpha_k k +\beta_k\\
\displaystyle
\rho\frac{\partial \varepsilon}{\partial t}  =
D(\varepsilon) + S_\varepsilon(k,\varepsilon)
+\varepsilon\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon)+\alpha_\varepsilon \varepsilon +\beta_\varepsilon
\end{array}\right.
\end{equation}

$D$ is the convection/diffusion operator.
$S_k$ (resp. $S_\varepsilon$) is the source term of $k$ (resp. $\varepsilon$).

\minititre{First step: explicit balance}

We solve the explicit balance:
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t} =
D(k^{(n)}) + S_k(k^{(n)},\varepsilon^{(n)})
+k^{(n)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n)})+\alpha_k k^{(n)} +\beta_k\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n)}) + S_\varepsilon(k^{(n)},\varepsilon^{(n)})
+\varepsilon^{(n)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n)})
+\alpha_\varepsilon \varepsilon^{(n)} +\beta_\varepsilon
\end{array}\right.
\end{equation}

(the term in $\Gamma$ is only taken into account in the case of forced injection)

\minititre{Second step: coupling of source terms}

The source terms are implicit in a coupled way:
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t} =
D(k^{(n)}) + S_k(k_{ts},\varepsilon_{ts})
+k^{(n)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n)})+\alpha_k k^{(n)} +\beta_k\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n)}) + S_\varepsilon(k_{ts},\varepsilon_{ts})
+\varepsilon^{(n)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n)})
+\alpha_\varepsilon \varepsilon^{(n)} +\beta_\varepsilon
\end{array}\right.
\end{equation}
so:
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t} =
\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}
+S_k(k_{ts},\varepsilon_{ts})-S_k(k^{(n)},\varepsilon^{(n)})\\
\displaystyle
\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}  =
\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
+S_\varepsilon(k_{ts},\varepsilon_{ts})-S_\varepsilon(k^{(n)},\varepsilon^{(n)})
\end{array}\right.
\end{equation}

The term in $\dive(\rho\vect{u})$ is not  implicit because it is linked to the
$D$ terms to ensure that the implication matrix
will be diagonally dominant. The term in $\Gamma$ and the user source
terms are not implicit either more, but they will be in the third step.

And we write (for $\varphi=k$ or $\varepsilon$)
\begin{equation}
S_\varphi(k_{ts},\varepsilon_{ts})-S_\varphi(k^{(n)},\varepsilon^{(n)})
=(k_{ts}-k^{(n)})
\left.\frac{\partial S_\varphi}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
+(\varepsilon_{ts}-\varepsilon^{(n)})
\left.\frac{\partial S_\varphi}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}
\end{equation}

So we finally solve the system $2\times 2$:
\begin{equation}
\left(\begin{array}{cc}
\displaystyle \frac{\rho^{(n)}}{\Delta t}
-\left.\frac{\partial S_k}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
&\displaystyle
-\left.\frac{\partial S_k}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}\\
\displaystyle
-\left.\frac{\partial S_\varepsilon}{\partial k}\right|_{k^{(n)},\varepsilon^{(n)}}
&\displaystyle
\displaystyle \frac{\rho^{(n)}}{\Delta t}
-\left.\frac{\partial S_\varepsilon}{\partial \varepsilon}\right|_{k^{(n)},\varepsilon^{(n)}}
\end{array}\right)
\left(\begin{array}{c}
(k_{ts}-k^{(n)})\\(\varepsilon_{ts}-\varepsilon^{(n)})
\end{array}\right)
=\left(\begin{array}{c}
\displaystyle\rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}\\
\displaystyle\rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
\end{array}\right)
\end{equation}

\vspace*{0.2cm}

\minititre{Third step: implicit convection/diffusion}

We solve the system:
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k^{(n+1)}-k^{(n)}}{\Delta t} =
D(k^{(n+1)}) + S_k(k_{ts},\varepsilon_{ts})
+k^{(n+1)}\dive(\rho\vect{u})+\Gamma(k_i-k^{(n+1)})\\
\multicolumn{1}{r}{+\alpha_k k^{(n+1)} +\beta_k}\\
\displaystyle
\rho^{(n)}\frac{\varepsilon^{(n+1)}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n+1)}) + S_\varepsilon(k_{ts},\varepsilon_{ts})
+\varepsilon^{(n+1)}\dive(\rho\vect{u})
+\Gamma(\varepsilon_i-\varepsilon^{(n+1)})\\
\multicolumn{1}{r}{+\alpha_\varepsilon \varepsilon^{(n+1)} +\beta_\varepsilon}
\end{array}\right.
\end{equation}
soit
\begin{equation}
\left\{\begin{array}{l}
\displaystyle
\rho^{(n)}\frac{k^{(n+1)}-k^{(n)}}{\Delta t} =
D(k^{(n+1)})-D(k^{(n)})+\rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}
+(k^{(n+1)}-k^{(n)})\dive(\rho\vect{u})\\
\multicolumn{1}{r}{-\Gamma(k^{(n+1)}-k^{(n)})+\alpha_k(k^{(n+1)}-k^{(n)})}\\
\displaystyle
\rho^{(n)}\frac{\varepsilon^{(n+1)}-\varepsilon^{(n)}}{\Delta t}  =
D(\varepsilon^{(n+1)})-D(\varepsilon^{(n)})
+\rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}
+(\varepsilon^{(n+1)}-\varepsilon^{(n)})\dive(\rho\vect{u})\\
\multicolumn{1}{r}{-\Gamma(\varepsilon^{(n+1)}-\varepsilon^{(n)})
+\alpha_\varepsilon(\varepsilon^{(n+1)}-\varepsilon^{(n)})}
\end{array}\right.
\end{equation}

The term in $\Gamma$ is not yet taken into account except in the case of
forced injection. The term in $\alpha$ is only taken into account if $\alpha$ is
negative, to avoid weakening the diagonal of the matrix that we are going to
invert.

\minititre{Coupling details}
During the coupling step, in order to favor stability and
realizability of the result, not all terms are taken into account
account. More precisely, we can write:

\begin{equation}
\left\{\begin{array}{l}
\displaystyle
S_k =
\rho C_\mu\frac{k^2}{\varepsilon}\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right)
-\frac{2}{3}\rho k \dive(\vect{u})
-\rho\varepsilon\\
\displaystyle
S_\varepsilon =
\rho C_{\varepsilon_1} C_\mu k\left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
-\frac{2}{3}C_{\varepsilon_1}\rho \varepsilon \dive(\vect{u})
-\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k}
\end{array}\right.
\end{equation}

Noting
$\displaystyle\tilde{\mathcal{P}}
= \left(\frac{\partial u_i}{\partial x_j} +
\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}
-\frac{2}{3}(\dive\vect{u})^2$\\
and
$\displaystyle\tilde{\mathcal{G}}
= -\frac{1}{\rho\sigma_t}
\frac{\partial\rho}{\partial x_i}g_i$

We therefore have in theory
\begin{equation}
\left\{\begin{array}{l}
\displaystyle \frac{\partial S_k}{\partial k}=
2\rho C_\mu\frac{k}{\varepsilon}\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right)
-\frac{2}{3}\rho \dive(\vect{u})\\
\displaystyle \frac{\partial S_k}{\partial \varepsilon}= -\rho\\
\displaystyle \frac{\partial S_\varepsilon}{\partial k}=
\rho C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
+\rho C_{\varepsilon_2}\frac{\varepsilon^2}{k^2}\\
\displaystyle \frac{\partial S_\varepsilon}{\partial \varepsilon}=
-\frac{2}{3}C_{\varepsilon_1}\rho \dive(\vect{u})
-2\rho C_{\varepsilon_2}\frac{\varepsilon}{k}
\end{array}\right.
\end{equation}

In practice, we will try to ensure $k_{ts}>0$ and $\varepsilon_{ts}>0$. In itself
based on a simplified calculation, as well as on feedback
of ESTET, we show that it is preferable not to take into account
certain terms. In the end, we realize the following coupling:

\begin{equation}
\left(\begin{array}{cc}
A_{11}&A_{12}\\
A_{21}&A_{22}
\end{array}\right)
\left(\begin{array}{c}
(k_{ts}-k^{(n)})\\(\varepsilon_{ts}-\varepsilon^{(n)})
\end{array}\right)
=\left(\begin{array}{c}
\displaystyle\frac{k_e-k^{(n)}}{\Delta t}\\
\displaystyle\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}
\end{array}\right)
\end{equation}
with
\begin{equation}
\left\{\begin{array}{l}
\displaystyle A_{11}=\frac{1}{\Delta t}
-2 C_\mu\frac{k^{(n)}}{\varepsilon^{(n)}}
\Min\left[\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right),0\right]
+\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]\\
\displaystyle A_{12}= 1\\
\displaystyle A_{21}=
- C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
- C_{\varepsilon_2}\left(\frac{\varepsilon^{(n)}}{k^{(n)}}\right)^2\\
\displaystyle A_{22}=\frac{1}{\Delta t}
+\frac{2}{3}C_{\varepsilon_1}\Max\left[\dive(\vect{u}),0\right]
+2 C_{\varepsilon_2}\frac{\varepsilon^{(n)}}{k^{(n)}}
\end{array}\right.
\end{equation}

(by definition of $C_{\varepsilon_3}$,
$\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}$
is always positive)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\etape{Calculation of the production term}
We call \fort{cs\_field\_gradient\_vector} to calculate the gradients of velocity.
In the end, we have \\
$\displaystyle \var{tinstk}=
2\left(\frac{\partial u}{\partial x}\right)^2+
2\left(\frac{\partial v}{\partial y}\right)^2+
2\left(\frac{\partial w}{\partial z}\right)^2+
\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2+
\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)^2+
\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2$\\
et\\
$\displaystyle \var{divu}=
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}
+\frac{\partial w}{\partial z}$

(the term $div(\vect{u})$ is not calculated by \fort{cs\_divergence}, for
correspond exactly to the trace of the strain tensor which is
calculated for production)

\etape{Reading user source terms}
Call \fort{cs\_user\_source\_terms} to load user source terms. They are
stored in the following tables:\\
$\var{W7}=\Omega\beta_k$\\
$\var{W8}=\Omega\beta_\varepsilon$\\
$\var{usimpk}=\Omega\alpha_k$\\
$\var{usimpe}=\Omega\alpha_\varepsilon$

Then we add the term in $(div\vect{u})^2$ \` to \var{TINSTK}. So we have \\
$\var{TINSTK}=\tilde{\mathcal{P}}$

\etape{Calculation of the gravity term}
Calculation only if $\var{igrake}=1$.\\
We call \fort{cs\_gradient\_scalar} for \var{rom}, with boundary conditions
$\var{coefa}=\var{romb}$ and \mbox{$\var{coefb}=\var{viscb}=0$}.\\
$\var{prdtur}=\sigma_t$ is set to 1 if there is no temperature scalar.

$\tilde{\mathcal{G}}$ is calculated and the source terms are updated:\\
$\var{tinstk}=\tilde{\mathcal{P}}+\tilde{\mathcal{G}}$\\
$\var{tinste}=\tilde{\mathcal{P}}+\Max\left[\tilde{\mathcal{G}},0\right]
=\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}$

If $\var{igrake}=0$, we simply have\\
$\var{tinstk}=\var{tinste}=\tilde{\mathcal{P}}$

\etape{Calculation of the mass accumulation term}
We store
$\displaystyle \var{W1}=\Omega\dive(\rho\vect{u})$
calculated by \fort{divmas} (to correspond to the convection terms of the
matrix).

\etape{Calculation of explicit source terms}
We assign the explicit source terms of $k$ and $\varepsilon$ for the
first step.\\
$\displaystyle\var{smbrk}=\Omega\left(\mu_t(\tilde{\mathcal{P}}+\tilde{\mathcal{G}})
-\frac{2}{3}\rho^{(n)} k^{(n)}\dive{\vect{u}}
-\rho^{(n)} \varepsilon^{(n)}\right)+\Omega k^{(n)}\dive(\rho\vect{u})$\\
$\displaystyle\var{smbre}=\Omega\frac{\varepsilon^{(n)}}{k^{(n)}}
\left(C_{\varepsilon_1}\left(
\mu_t(\tilde{\mathcal{P}}+(1-C_{\varepsilon_3})\tilde{\mathcal{G}})
-\frac{2}{3}\rho^{(n)} k^{(n)}\dive{\vect{u}}\right)
-C_{\varepsilon_2}\rho^{(n)}\varepsilon^{(n)}\right)
+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})$

then $\var{smbrk}=\Omega S_k^{(n)}+\Omega k^{(n)}\dive(\rho\vect{u})$
and $\var{smbre}=\Omega S_\varepsilon^{(n)}+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})$.


\etape{Calculation of user source terms}
We add the explicit user source terms to \var{smbrk} and
\var{smber}:\\
$\var{smbrk}=\Omega S_k^{(n)}+\Omega k^{(n)}\dive(\rho\vect{u})+\Omega\alpha_k k^{(n)} +\Omega\beta_k$\\
$\var{smbre}=\Omega S_\varepsilon^{(n)}+\Omega\varepsilon^{(n)}\dive(\rho\vect{u})
+\Omega\alpha_\varepsilon \varepsilon^{(n)} +\Omega\beta_\varepsilon$

Arrays \var{w7} and \var{w8} are freed, \var{usimpk} and \var{usimpe} are
retained for use in the third step of resolution.

\etape{Computation of explicit convection/diffusion terms}
\fort{cs\_balance\_scalar} est called twice, for $k$ et for $\varepsilon$, so
as to add to \var{smbrk} and \var{smbre} the explicit convection/diffusion terms
explicites $D(k^{(n)})$ et $D(\varepsilon^{(n)})$. These termes are first
stored in \var{w7} and \var{w8}, to be reused in the third resolution stage.

\etape{Mass source terms}
In the case of a forced injection of matter, we call \fort{ cs\_mass\_source\_terms}
twice to add the terms in
$\Omega \Gamma (k_i-k^{(n)})$ and
$\Omega \Gamma (\varepsilon_i-\varepsilon^{(n)})$ to \var{smbrk} and
\var{smber}. The implicit part ($\Omega\Gamma$) is stored in the
variables \var{w2} and \var{w3}, which will be used during the third
step (the two variables are indeed necessary, in case we have a
forced injection on $k$ and not on $\varepsilon$, for example).

\etape{End of the first step}
This completes the first step. We have\\
$\displaystyle \var{smbrk}=\Omega \rho^{(n)}\frac{k_e-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{smbre}=\Omega \rho^{(n)}\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}$

\etape{Coupling step}
(only if $\var{ikecou}=1$)

We renormalize \var{smbrk} and \var{smbre} which become the rhight hand side terms of the
coupling system.\\
$\displaystyle \var{smbrk}=\frac{1}{\Omega\rho^{(n)}}\var{smbrk}
=\frac{k_e-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{smbre}=\frac{1}{\Omega\rho^{(n)}}\var{smbre}
=\frac{\varepsilon_e-\varepsilon^{(n)}}{\Delta t}$\\
and $\displaystyle \var{divp23}=\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]$.

We set the coupling matrix\\
$\displaystyle \var{A11}=\frac{1}{\Delta t}
-2 C_\mu\frac{k^{(n)}}{\varepsilon^{(n)}}
\Min\left[\left(\tilde{\mathcal{P}}+\tilde{\mathcal{G}}\right),0\right]
+\frac{2}{3}\Max\left[\dive(\vect{u}),0\right]$\\
$\displaystyle \var{A12}= 1$\\
$\displaystyle \var{A21}=
- C_{\varepsilon_1} C_\mu \left(\tilde{\mathcal{P}}
+(1-C_{\varepsilon_3})\tilde{\mathcal{G}}\right)
- C_{\varepsilon_2}\left(\frac{\varepsilon^{(n)}}{k^{(n)}}\right)^2$\\
$\displaystyle \var{A22}=\frac{1}{\Delta t}
+\frac{2}{3}C_{\varepsilon_1}\Max\left[\dive(\vect{u}),0\right]
+2 C_{\varepsilon_2}\frac{\varepsilon^{(n)}}{k^{(n)}}$

We invert the system $2\times 2$, and we get:\\
$\displaystyle \var{deltk}=k_{ts}-k^{(n)}$\\
$\displaystyle \var{delte}=\varepsilon_{ts}-\varepsilon^{(n)}$

\etape{End of the second step}
We update the variables \var{smbrk} and \var{smbre}.\\
$\displaystyle \var{smbrk}=\Omega \rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}$\\
$\displaystyle \var{smbre}=
\Omega \rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}$
If we do not couple ($\var{ikecou}=0$), these two variables keep
the same value as at the end of the first step.

\etape{Calculation of implicit terms}
We remove from \var{SMBRK} and \var{SMBRE} the part in convection diffusion at
time $n$, which was stored in \var{W7} and \var{W8}.\\
$\displaystyle \var{SMBRK}=\Omega \rho^{(n)}\frac{k_{ts}-k^{(n)}}{\Delta t}
-\Omega D(k^{(n)})$\\
$\displaystyle \var{SMBRE}=
\Omega \rho^{(n)}\frac{\varepsilon_{ts}-\varepsilon^{(n)}}{\Delta t}
-\Omega D(\varepsilon^{(n)})$

We calculate the implicit terms, excluding convection/diffusion, which correspond
diagonal matrix.\\

$\displaystyle \var{tinstk}=\frac{\Omega \rho^{(n)}}{\Delta t}
-\Omega\dive(\rho\vect{u})+\Omega\Gamma+\Omega\Max[-\alpha_k,0]$\\
$\displaystyle \var{tinste}=\frac{\Omega \rho^{(n)}}{\Delta t}
-\Omega\dive(\rho\vect{u})+\Omega\Gamma+\Omega\Max[-\alpha_\varepsilon,0]$\\
($\Gamma$ is only taken into account in forced injection, i.e. it
is necessarily positive and does not risk weakening the diagonal of the matrix).

\etape{Final Resolution}
We then pass twice in the function \fort{cs\_equation\_iterative\_solve\_scalar}, for $k$ and $\varepsilon$,
to solve equations of the type:

$\var{tinst}\times(\varphi^{(n+1)}-\varphi^{(n)}) = D(\varphi^{(n+1)})+\var{smbr}$.

\etape{Final clipping}
We finally pass in the function \fort{cs\_turbulence\_ke\_clip} to clip
clip $k^{(n+1)}$ et $\varepsilon^{(n+1)}$ if necessary.
